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Here we give detailed derivations and provide additional examples to the main paper In 
particular, we discuss the scaling behavior of observables like correlation functions and density of 
excitations. We also analyze effects of nonintegrability of the Bose-Hubbard model on the long-time 
dynamics of the correlation functions. In addition we explicitly consider several interacting models, 
where we are able to analyze slow dynamics and classify it according to the regimes suggested in 
the main paper. 
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I. INTRODUCTION AND DISCUSSION 

The aim of this Supplementary Material is twofold. 
First, in Sections II- V we provide more details to the 
main text dj. Thus in Sec. [II] we describe the approach 
to the slow dynamics in our model based on the Fermi 
Golden Rule in the ramping rate. In Sec. IIIII we give 
details of the time evolution of the wave function of the 
harmonic model if the system is initially prepared in the 
ground state. Then in Sec. IIVI we generalize this deriva- 
tion to the evolution of the density matrix assuming 
the initial thermal distribution. Dynamics of the Bose- 
Hubbard model is addressed in Sec.lVl There we give the 
details of our numerical approach based on the semiclassi- 
cal approximation and the leading quantum corrections. 
We also discuss the consequences of non-integrability of 
the Bose-Hubbard model on the time evolution of the 
correlation functions. In particular, we show that at very 
long times the non-equilibrium state created during the 
ramp relaxes to the thermal equilibrium. 

In the main paper [l[ we exclusively concentrated on 
finding energy added to the system during the ramp. In 
this Supplementary Material we will also consider vari- 
ous other quantities like density of excitations, correla- 
tion functions etc. Sometimes these quantities are eas- 
ier to measure experimentally and if the effects of non- 
integrability are weak or absent then they are good ob- 
servables to work with. 

The second aim of this Supplementary Material is to 
consider application of our findings to several physical 
problems in more details. In Sec. [VI] we discuss poten- 
tial relevance of our findings to the quantum information 
and to some problems in inflationary cosmology. Then in 
Sec. IVIll discuss the dynamics of one dimensional bosons 
in Tonks- Girardeau regime 0. In particular, we show 
that if the trapping potential for atoms is slowly reduced 



to zero then the heating induced in the system is de- 
scribed by the non-analytic B) regime according to our 
classification. We show that these results can also be 
applied to the Calogero-Sutherland model @, describing 
one dimensional fermions with long range interactions, 
in the harmonic trap. And finally in Sec. I Villi we briefly 
describe time evolution of the quantum Dicke model 0] , 
which serves as a prototype for the laser, as well as mim- 
ics a coherent atomic cloud in the cavity QED. There in 
a particular regime we find that the number of generated 
photons is described by the non-adiabatic regime C). Of 
course this list of possible applications of our findings is 
not complete, yet it is quite illuminating. We note that 
some of these models are integrable, some are not. Yet 
we find general good agreement of the slow dynamics in 
all of these models with our classification scheme. 



II. FERMI GOLDEN RULE ANALYSIS OF THE 
SLOW DYNAMICS. 



One way to find the density of excitations n cx and 
the energy produced during a slow increase is to use the 
Fermi Golden in the ramp speed Q . We remind that we 
consider the Hamiltonian (5) from the main text: 
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where we choose K q = k + Xq and k = kq + St linearly 
change in time. Since n ex should be small at small 5 one 
can expect that the perturbation theory in S gives a good 
estimate of n ex (6). Then the density of excitations n ex 
can be expressed as follows @ : 
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where \m) K denotes a general excited many-body state 
with the energy Two m (K,), (m\d/ cLk\0) K is the matrix el- 
ement of the derivative with respect to k between the 
states |m) and |0) at a given value of n. We would like to 
emphasize that Eq. ^ is valid only if the ground state 
evolution does not acquire an additional Berry phase. In 
the situation where the Berry phase is nonzero it should 
be subtracted from the argument of the exponent in this 
equation. It is straightforward to check that with the 
Hamiltonian ([1]) the only non-vanishing matrix element 
of d/dn corresponds to the excitation of two quasiparti- 
cles with opposite momenta: (q, — q|<9 K |0) = 1/(4v / 2k ( j). 
Using that u) m — ujq = 2q^jp s n qi where the factor of 
two comes from the fact that we have two excited quasi- 
particles, we find: 
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where ko, q = kq + Xq 2 . This expression gives different 
asymptotics in t he tw o opposite limits. 

(i) If S 3> Kq \/ Ps/X, which is the case if one starts from 
the weakly interacting regime kq — > 0, then 
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where Ad is a numerical constant. It is easy to check 
that if d > d* = 8 the exponent of 8 saturates at 2 and 
does not depend on the dimensionality. Expression (|4|) 
suggests that in this particular situation the nonanalytic 
regime B is realized in all physical dimensions. In one 
dimension it is particularly hard to reach the adiabatic 
regime since n ox scales only as 5 1 / 4 . 

We note that the scaling in Eq. is consistent with 
the one obtained in Ref. [5( for the crossing of the second 
order phase transition: n cx cx S dv / i zv + x ) ^ where v is the 
critical exponent characterizing divergence of the corre- 
lation length. In our case there is a diverging healing 
length £ ~ \/X/k instead of the correlation length (see 
Ref. @ for details) so that v = 1/2 and given that z = 2 
in the noninteracting regime one immediately recovers 
that v/(zv+ 1) = 1/4. 

(ii) In the opposite limit, where the initial value of k is 
large 5 -C kq^J p s /\ the situation becomes more diverse. 
Thus for dimensions d < 2 Eq. |(3} yields 
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On the other hand for d > 2 the exponent saturates and 
we have 
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In two dimensions there is an additional logarithmic cor- 
rection to the scaling ([5]). We see that in this situation 



the critical dimension above which the analytic regime 
holds is d* = 2. 

The present analysis can be generalized to other situ- 
ations. For example, in the case of ferromagnets kq = 
and then one can tune A. Then one finds that n cx oc 
and the critical dimension is d* = 4. We comment that 
one can also consider other scenarios of varying k with 
time. For example, if k oc (5t) r then it is easy to see that 
n GX oc (5 <ir / 2 ( r + 1 ) . As r increases the scaling of the density 
of excitations interpolates from <5 d//4 to 5 d ^ 2 and changes 
d* from eight to four consistent with a recent prediction 
of Ref. 0- 

This perturbative analysis shows the existence of A 
(B) regimes for dimensions above (below) some critical 
value d* . However, it misses the existence of the non- 
adiabatic C regime. To justify the validity of the appli- 
cation of the Fermi golden rule one has to require that 
the probability of excitation of each momentum mode 
is small. This requirement breaks down at low energies 
as can be readily seen from Eq. ([3]). In the case when 
the excitations have Fermionic character, which is e.g. 
the case for crossing the critical point in the transverse 
field Ising model or the XXZ chain [f| , the mistake of the 
perturbative treatment is a simple factor of the order of 
one (see Refs. 0, H [lfl E2). The Goldstone modes de- 
scribed by the Hamiltonian (fl]) are harmonic oscillators 
and thus behave as bosons. Bosons unlike fcrmions have 
a bunching tendency, i. e. transition probabilities can 
be significantly enhanced compared to the golden rule 
prediction. 

For the energy density in the system one can derive a 
similar expression to Eq. ^ : 
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From here, for example, for the initially noninteracting 
case kq = one recovers Eq. (6) of the main text: 
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Similarly one can reproduce the correct scaling for finite 
kq mentioned in the main text. 



III. EVOLUTION OF THE WAVE FUNCTION 
AT ZERO INITIAL TEMPERATURE. 

The harmonic theory described by the Hamiltonian ([1]) 
can be in principle analyzed for arbitrary functional de- 
pendence of coupling K q (t) . This is a consequence of the 
fact that the quantum harmonic oscillator problem can 
be solved for arbitrary functional dependence of its fre- 
quency (and mass) on time. The resulting Riccati-type 
equation can be analytically solved in some situations. In 
particular, this is the case for the linear time dependence, 
which we analyze here in more detail. 
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As we described in the main text the initial ground 
state wave function is given by 
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where ao, q — l/(2q)^/ Ko,q/ Ps- If K changes with time, 
o q acquires time dependence: 
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This equation can be simplified by first changing indepen- 
dent variable t to K q (t) and then by a simple rescaling: 



52/3 

K- 



51/3 



51/4 



1/8 



A 3 / 8 



Under these transformations we also have k q = k 4 
Then one can check that Eq. (JTUJ) is equivalent to 

.da „ 



dk n 
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This Riccati equation can be explicitly solved in terms of 
Airy functions Ai and Bi: 



. Bi(-k q ) + a q Ai'(-k q ) 
Bi(-K g ) + a q Ai(-k q ) 



(13) 



where a q is an integration constant, which is determined 
from the initial conditions. In the limit k q - 
unimportant fast oscillating terms we find 
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Note that the real part of l/a q determines \ifj\ 2 and thus 
the probability distribution of the corresponding Fourier 
component of the phase (f> q (see Eq. (|9])). The fact that 
1/fq ~~ * as K q — ► 00 should not be surprising. Indeed 
the width of the ground state wave function in scaled 
variables is 



(15) 



The probability of excitations in the system is determined 
by the ratio of a q and tr eq , which takes a well defined limit 
at k — > 00 . Introducing a q s = l/SJf^cj" 1 ) we find 
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The initial condition determining a is: 

.Bi'(-Ko, 9 ) + a g Ai'(-Ko,g) 



Bi(-Ko g ) + a q Ai(-k 0:q ) 



(16) 



(17) 



This equation can be inverted to give 

Jk . q Bi(-K , q ) - iBi'(-k q ) 
y/Ro, q Ai(-KQ, q) - iAi'(-«o, q) 

In the limit Ko i9 <C 1 this equation yields: 
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In the opposite limit Kq 3> 1 one finds 
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We note that in this limit Eq. (|21[) gives the result iden- 
tical to what one would get using Fermi Golden rule ap- 
proach described in the previous section: 
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We remind that n q and a q are related according to 
Eq. (13) of the main text: 
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so that Eqs. (|2"Tj) and (f2"2"|) indeed agree for high energy 
modes. 

The number of excitations studied above is not neces- 
sarily an observable quantity. Instead one can look, for 
example, into the behavior of the correlation functions, 
which are closely related to the population of different 
modes: 
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If the initial state is noninteracting: kq = then ac- 



cording to Eqs. (f2"0"l) . (f]~5|) . and pip we have <r* ff oc q 
at small q. Therefore in one and two dimensions 
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(25) 



In one dimension this integral decays faster than expo- 
nential indicating that the system is overheated, i.e. the 
behavior of the correlation functions is worse than at fi- 
nite temperature. In two dimensions the correlation func- 
tions decay as exp[— CS^-^x 1 ^ 3 ], which is again a very 
unusual behavior. Note that in one and two dimensions 
the asymptotic behavior of correlation functions (|25D is 
valid only at long distances x > £<j with ~ 1/5 1 / 4 



4 



and ^2D ~ 1/5. In dimensions d > 7/3 the excitations 
in the system do not destroy the long-range order in the 
system but reduce the superfluid density 
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If one starts in the interacting regime: kq 3> y <5/no 
then one finds that <7 c ff(<7) ^ at small g and thus 

the correlation functions are singular only in one dimen- 
sion: 
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In this case the correlation length diverges as £m ~ 1/8. 
Above one dimension the long range order survives and 
the long distance behavior of correlation functions is: 
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for d < 3. Above three dimensions the power of S in the 
expression above saturates at two. 

We would like to stress that the steady state non- 
equilibrium distribution of quasi-particles and as a con- 
sequence noneqiulibrium correlation functions, which we 
obtained above are only possible in strictly noninteract- 
ing model. Indeed addition of small nonlinear terms into 
the Hamiltonian ([1]) can lead to redistribution of excita- 
tions among different states and eventual thermalization. 
We already highlighted that this is indeed the case in the 
main text and will return to this issue again in Sec. [V] 
We note that these possible thermalization processes do 
not affect the total energy, which is conserved in an iso- 
lated system. 



IV. EVOLUTION OF THE DENSITY MATRIX 
AT FINITE INITIAL TEMPERATURE. 

We choose to represent the density matrix corre- 
sponding to the initial thermal state in the Wigner 
form jl2l [l3T | . For the harmonic system described by the 
Hamiltonian (5) one can show that this density matrix 
factorizes into the product of Gaussian functions: 
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In the noninteracting problem the time evolution of the 
fields 4>q and Iln is described by the classical equations 
of motion [13l \i 



subject to the initial conditions 

<j) q {t = 0) = 0o, 9 , 4>g(t = 0) = Ko^IIo^. (32) 

Here (j)Q, q and Ho, q are randomly distributed accord- 
ing to Eq. (p9|) . The other important feature of Gaus- 
sian ensembles is that in the absence of interactions the 
Wigner distribution (|29| always preserves its Gaussian 
form. Therefore finding (<f> q (t)) and (U^(t)) is sufficient 
to fix the whole distribution function at arbitrary time. 
Alternatively one can directly solve the Liouville equa- 
tion for the density matrix in the Wigner form fl3j and 
come to the same conclusion. 
A general solution of Eq. ([3~l]) is: 



<f> q (kq) = CiAi'(-K q ) + C 2 Bi'(-K g ), 



(33) 



where as in Appendix [Hi] we changed the variables from t 
to k q . The integration constants C\ and C2 can be found 
from the initial conditions: 



Ci 



— — -7-^Bi'(-K , g ) - 7r0o, g Bi(-Ko, 9 ), (34) 

K 0, <j " ,K 0, q 

C 2 = 7T0O, gAi(-K , q) ~ Ai'(-K , q ). (35) 
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From these expressions it is easy to find the asymptotical 
behavior of (0 2 ) at large k and thus find the width of the 



distribution <j^ ff : 



-!q = 77 ?S — [^0. ? Bi 2 (-K ,g) + K0,gAi 2 (-K . g ) 

u q i y^o, q 

+ (Bi'(-« , 9 )) 2 + (Ai'(-^ g )) 2 ]. (36) 

One can verify that apart from the factor of r q , which 
approaches unity at T — > 0, the expression above coin- 
cides with the zero temperature results (see Eqs. (116[) and 
(HE])), i.e. 
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This result immediately implies that the number of the 
additional excitations created during the ramp at finite 
temperature can be obtained from the zero temperature 
result by multiplication by r q : 



1 



r q n q \ T=Q + ~(r q - 1) 



(38) 



Integrating n q over momenta we find that the total 
density of excitations for kq = scales at finite temper- 
ature in all three spatial dimensions as 



n cx oc TL 10 / 3 - d V6. 



(39) 
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(31) 



So in terms of the excitation density the system is always 
in the regime C). The energy density shows less divergent 
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behavior and the regime C) is realized only in one and 
two spatial dimensions while in the three dimensional 
case dynamics belongs to the nonanalytic B) as described 
in the main text 

If the initial compressibility kq is finite than we find 
that in one dimension the density of excitations still di- 
verges with the system size: 



oc Ty 5L, 



(40) 
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but it is finite in two and three dimensions rt cx oc TS 
(as before the exponent of S saturates at two for d > 3). 
The total energy converges in all three dimensions and it 
behaves as £f oc T6 d for d < 2 and £f cx TS 2 for d > 2. 

As in the previous section one can compute correla- 
tion functions. Note that because for initially noninter- 
acting regime a q diverges as 1/g 13 / 3 the sum in Eq. ([24]) 
is infrared divergent in one and two dimensions even at 
qx < 1. This results in a very unusual behavior of the 
correlation functions. 
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In three dimensions we have 

(e^-^°») ~ exp [-CTVSx 4 / 3 ] . (42) 

We comment again that this unconventional behavior of 
the correlation functions can exist as long as the relax- 
ation processes in the system are negligible. We will get 
back to this issue in the next section. 



V. QUANTUM DYNAMICS OF A 
BOSE-HUBBARD SYSTEM: EXPANSION IN 
QUANTUM FLUCTUATIONS. 

Here we will describe in some detail how to simulate 
slow dynamics of the system described by the Hubbard 
model (11) of the main text using the semiclassical ap- 
proach [14J. For completeness we will write the Bosc- 
Huabbard Hamiltonian again: 



n bh = -j£>k + a]ai) + 

(ij) 3 

(43) 

Here a,j and al are the bosonic annihilation and creation 
operators, J represents the tunneling matrix element and 
U is the interactions strength. The sum in the first term 
is taken over the nearest neighbor pairs. 

Specifically we will use expansion of the time evolution 
of the system in the small quantum parameter U/Juq. 
Note that when this parameter is close to one, the ground 
state of the system undergoes the superfluid-insulator 
transition driven by quantum fluctuations [§[ . Conversely 
when U/JriQ <C 1 quantum fluctuations are negligible 
and the system is in the superfluid regime. From this 
one can conclude that this ratio plays the role of the 



Planck's constant in this problem (see Ref. [15j for more 
details). Here we are interested in evolution precisely 
in the regime where the system is far from the insulating 
phase and the harmonic approximation is accurate so the 
expansion in this ratio is justified. 

We note for those more familiar with the Kcldysh tech- 
nique [16j that our approach treats all classical vertexes 
exactly and expands the evolution in number of quan- 
tum vertexes. In the leading order in this parameter 
one obtains the so called truncated Wigner approxima- 
tion (TWA) [12, [n}, where the classical fields ip*- and 

ipj corresponding to the operators aj and ctj satisfy the 
time dependent Gross-Pitaevskii equations of motion. In 
the next order the classical fields are subject to a single 
quantum jump during the evolution. We find that while 
TWA approximation is adequate at finite temperatures, 
in the zero temperature limit one has to go beyond and 
add the next correction. This finding agrees with a gen- 
eral statement that the semiclassical approximation can 
break down at long times [3, [l|| . 

In the classical limit the bosonic fields ipj and ipj satisfy 
the discrete Gross-Pitaevskii equations: 



.dip 



^± = -jj2A + u(t)\ip^ 



(44) 
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Here the sum in the first term is taken over the nearest 
neighbors of the site j. In the leading order in quan- 
tum fluctuations, which corresponds to the semiclassical 
or truncated Wigner approximation (TWA) , the fields ipj 
and ipj are subject to random initial conditions, which 
are distributed according to the Wigner transform of the 
initial density matrix W(ip* ,ipj). The expectation value 

of an arbitrary observable fi(aj, a,j) is given by the aver- 
age of the corresponding Weyl symbol (fully symmetrized 
form of the operator) fl c \(ipj, ipj) on the solutions of the 
Gross-Pitaevskii equations: 



(fi(t)) = J D^DiP.W^^j^cM^mm- (45) 

Since the initial system is noninteracting, it is straight- 
forward to find the Wigner transform of the density ma- 
trix at finite temperature T. It is more convenient to 
write it in the Fourier space 



W(tPt,ip k ) = Zjjexp 



-2|^,| 2 tanh 



ep(ff) - M 
2T 



(46) 

where ipk is the discrete Fourier transform of ipj , Z is the 
normalization constant, eo(q) — —Jj2j ^ i s the excita- 
tion energy of the Bose-Hubbard Hamiltonian (|43[) in the 
absence of interactions and the summation is taken over 
nearest neighbors of site at the origin, is the chemical 
potential which enforces mean number of particles per 
site no. We note that in large systems we consider here, 
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there is no difference in time evolution between grand 
canonical and canonical ensembles [Tl |. 

We find that the semiclassical approximation (j45| gives 
very accurate results in most of our simulations described 
in this paper. However, at zero temperature case it 
breaks down for very slow ramps and we had to include 
the next quantum correction to the TWA. The latter 



manifests itself in the form of a single infinitesimal quan- 
tum jump during the evolution: 

tp i {t')^i) i {t') + e l +ie 2 . (47) 

The quantum correction is the evaluated as a nonlinear 
response of f2 c i to such a jump [i"4| : 



m))i = - / d^d^w^i^y 



1G 



5hl>iV)-S--W>iV)-S- 
oei de 2 





" d 2 


d 2 - 






h de 2 _ 



n cl (^*(t),^(t),e 1 ,e 2 ). 
(481 



Numerically both the leading term (£l(t))o and the 
next correction (f2(t))i are evaluated using Monte-Carlo 
integration schemes. The third order derivatives in 
Eq. (|48|) are found using finite differences, e. g. 



d 3 n(ei) »(2ei) - fl(-2ei) - 2»(ei) + 2»( £l ) 

del ~ 
d 3 n(e 1 ,e 2 ) 



2e\ 



(49) 



de\de\ 



... ll O(ei,e2) + 0(e 1 ,-e 2 ) (50) 
Zeie 2 



-n(-ei, e 2 ) - fi(-ei, -e 2 ) - 2fi(ei,0) + 2fi(- ei) 0) 



It is easy to convince oneself that in order to evaluate 
these finite differences one has to simultaneously solve 
thirteen Gross-Pitaevskii equations, one for £1 = 62 = 
and the others for various combinations of £1,62 = 
0, ±e, ±2e. While solving thirteen Gross-Pitaevskii equa- 
tions is certainly more time consuming task than solving 
one equation, it is still tremendously more advantageous 
than dealing with the exact quantum problem. To il- 
lustrate the importance of quantum correction at zero 
temperature we show comparison of dependence AS (6) 
at zero temperature with and without this correction (see 
Fig. [T]). The semiclassical approximation gives spurious 
saturation (and even increase) of the heating induced in 
the system as S — ► 0. At the same time adding the first 
correction removes this unphysical behavior and extends 
the validity of the numerical results to slower rates. 

It is interesting that at finite (even very small) temper- 
atures the domain of validity of the semiclassical (TWA) 
approximation tremendously increases. Indeed the de- 
pendence of AS on 6 does not show any spurious behav- 
ior down to the slowest rate we were able to analyze (see 
Figs. 1-3 of the main text). This result is perhaps in- 
tuitively clear: we expect that quantum corrections play 
smaller role at higher temperatures. Nevertheless it is 
still quite surprising that even very small temperature 
T = 0.02, corresponding to only 1% of the band width 
2 J has such a strong effect on the validity of the semi- 
classical making it virtually exact. 

The next very important issue we would like to address 
here is whether the Bose-Hubbard model indeed leads 




FIG. 1: Dependence of the energy density A£ on the S at 
zero temperature with and without the quantum correction 
(|48[) . For the details of the calculation and the parameters of 
the problem see Fig. (1) of the main text. Obviously at small 
values of 5 the TWA breaks down and one has to include the 
correction 



to eventual thcrmalization and how it affects observables 
other than energy. We emphasize that the Bose-Hubbard 
model is not integrable in all spatial dimensions and 
thus thermalization is expected. However, at low ener- 
gies the excitations of this model are weakly interacting 
long- wavelength phonons. Thus on general grounds one 
can expect that the relaxation times of these phonons 
are very long. We also point that neglecting relaxation, 
during our process we primarily populate low energy ex- 
citations, which generically have longer life times than 
the high energy excitations. Thus we expect that the 
relaxation times in our case will be even longer than in 
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equilibrium. 

To analyze thermalization in the system we will con- 
centrate on the behavior of the correlation functions. Ac- 
cording to the noninteracting theory the long distance be- 
havior of these correlation functions is given by Eq. (|25p 
if the initial temperature is zero and by Eq. (j4"Tj) at a 
finite T. In Fig. [2] we plot correlation functions (a^ao) 



Correlation Functions 



Numerical 
Analytical 




5 10 15 20 

L/tc sin(7ry/L) 



25 



FIG. 2: Correlation function (ajdi+j) of a one-dimensional 
bosonic Hubbard model as a function of scaled distance 
L/-K sin(7Tj/L) at t = 3.2/5. The parameters of the model 
are identical to those in Fig. (4) of the main text: S = 0.1, 
L = 256, T = 0.02. The two lines represent the numer- 
ical data and the analytical result evaluated according to 
Eqs. U2U, (0, G3, and fl37}. 

at t = 3.2/(5 evaluated numerically (solid black line) and 
analytically according to Eqs. f2"3|). (fig]), (fig)) , and (jg?)) . 
Because we are dealing with a discrete system we need to 
change in all expressions q — > 2sin(q n /2) — 2sin(7rn/L), 
where n is an integer. The sum in Eq. (|24|) is taken over 
n = 1, . . . , L- 1. We use 6 = 0.1, L = 256, and T = 0.02 
- the same parameters as in Fig. (4) of the main text. 
The time is chosen such that the interaction U is almost 
saturated at U$ and yet the system did not have time to 
relax to the ground state. The agreement between the 
two curves is quite good, especially given the crudeness of 
the analytic approximation at this value of 5, where the 
heating is significant and the harmonic approximation to 
the Hubbard model is not expected to be very accurate. 
The small deviation between the two curves can be also 
due to partial relaxation of the system by the observation 
time. We emphasize that there are no fitting parameters 
involved in this comparison. We expect that the agree- 
ment between analytic and numerical results should be 
even better for smaller values of 5. 

As we mentioned above the system should eventually 
thermalize and the correlation functions should assume 
the equilibrium form. And indeed it happens as it is 
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FIG. 3: Correlation function (a]ai+j) of a one-dimensional 
bosonic Hubbard model as a function of scaled distance 
L/tt sin(7rjf/L) at different moments of time. Here 5 = 0.1, 
L = 256, T = 0.02. This figure duplicates Fig. 4of the main 
text. 



shown in Fig. 4 of the main text [l|, which we repeat 
here for completeness (Fig. [3]). The shape of the corre- 
lation function clearly evolves in time and approaches a 
steady state, which is very close to the thermal equilib- 
rium. We point out again that the thermal distribution 
is obtained for the noninteracting model with the tem- 
perature extracted from the total energy of the system 
and thus there are no fitting parameters involved. Obvi- 
ously the short distance part of the correlation functions 
thcrmalizcs faster than the its long distance tail. This 
observation is consistent with general expectations that 
thermalization times for short wavelength excitations are 
shorter. We emphasize that the overall relaxation time 
is very long of the order of 10 2 — 10 3 , while the natural 
time scales in the problem, like inverse Josephson fre- 
quency or the inverse frequency associated with steady 
state temperature T ~ 4.8, are much shorter indicating 
that the phonon relaxation times are very long. If one 
goes to smaller values of 5 then the thermalization time 
dramatically increases and the non-equilibrium shape of 
the correlation function can be observed for long times. 



VI. APPLICATION TO COSMOLOGY AND TO 
ADIABATIC QUANTUM COMPUTATION 

A. Small-roll approximation in cosmology. Particle 
creation. 

The problem of evolution of Universe is inherently adi- 
abatic in nature. The inflationary cosmology is essen- 
tially defined by the following equations for the scalar 
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field (f> which lives in a potential V(4>). 

+ 3ff0 + — =0, H 2 = — ( -4? + V) , (51) 

It is important that the potential V = V((f>) has a flat 
(constant) part. This allows us to introduce the so- 
called "slow-roll" approximation which can be defined 
as 2 < V(cj)), \4>\ < \H<f>\, and \4>\ < \8V/d<j>\. These 
approximations imply that the dynamics of the field <j) 
is slow, adiabatic. The whole scenario of the inflationary 
cosmology is based on this assumption (see e.g. Ref. (l9|). 
However, as we showed in this paper initial quantum or 
thermal fluctuations can be enhanced if the Universe is 
in the non-adiabatic regime. 

Other potential applications of our findings in cosmol- 
ogy include the problems of particle creation in the ex- 
panding Universe [2(J, the problem which brought a lot 
of attention in the literature (see e.g. Ref. [2l| and refer- 
ences therein). There, the quantum Hamiltonian of the 
fluctuating scalar field has essentially the same form as 
we used in our manuscript: 

H = \Y. l a 'AI 2 + ( fc2 - ^f) d >) (52) 

Here the conformal time r\ is rescaled by the scaling fac- 
tor a(t) which is responsible for the expansion in the 
Friedmann-Robertson- Walker metric. The specific ef- 
fects, which follow from application of our formalism 
to these problems remain to be investigated. We note 
that recently established connections between expansion 
of the Universe and of the Bose-condensate from a time- 
dependent trap 22] can be used to experimentally inves- 
tigate the effects of non-adiabaticity in the Universe. 

As it follows from our work, possible non-adiabatic ef- 
fects can be quite significant. In particular, the physics of 
cosmic microwave background radiation (CMB) includ- 
ing the prediction of the temperature is always described 
using the assumptions of adiabaticity of the expansion of 
the Universe. We hope these issues will be addressed in 
future by specialists working in cosmology. 

B. Adiabatic quantum computation. 

The concept of adiabatic quantum computation was 
originally proposed in Ref. [23[ as a method of solving 
combinatorial optimization problems. In this approach 
one starts with a quantum Hamiltonian for which the 
ground state can be easily constructed. Then the Hamil- 
tonian is adiabatically changed into another one, whose 
ground state encodes the solution of the problem. The 
use of the adiabatic theorem guarantees that the sys- 
tem will remain in the instantaneous ground state if the 
variation of the Hamiltonian is sufficiently slow. There 
has been a grown interest in using adiabatic quantum 
computation as an architecture for experimental quan- 
tum computation schemes. For practical applications of 



the adiabatic quantum computation it is very important 
whether this scheme has inherent fault tolerance. Un- 
derstanding this issue inspired interest to fundamental 
questions of the general applicability of the adiabatic the- 
orem 24]. Indeed it was argued that there might be an 
inconsistency or insufficiency of conditions in applicabil- 
ity of the adiabatic theorem and for some specific physical 
systems. These results further motivated large amount 
of works [H, [H, [13, [H, HI, IU HI examining the appli- 
cability of adiabaticity for variety of systems, including 
those which are envisioned for quantum computations 
and cavity QED. Usually in these works fidelity, i.e. the 
overlap of the wave function with the ground state, is 
used as a measure of non-adiabaticity. It is possible that 
real computational schemes can tolerate small number 
of excitations in the system. Our analysis suggests that 
even this weaker requirement of adiabaticity can be hard 
to achieve if regimes B) or C) are realized. 

VII. EVOLUTION OF AN INTERACTING ID 
GAS IN A TIME-DEPENDENT TRAP: 
APPLICATION TO THE TONKS GAS AND THE 
CALOGERO-SUTHERLAND MODEL. 

According to the results of Ref. [HJ the Tonks gas 
in an arbitrary time dependent parabolic trap V ex t = 
mu> 2 (t)x 2 /2 can be described exactly using the scal- 
ing approach. The evolution of the wave function of 
TG gas of N particles is given by the wave function 
®tg(xi, ■ ■ ■ , xn', 0) of the gas at initial time t = 

$tg(%i, ■■■,XN]t) = -£WJ5®TG{xi/b, ■ ■ .,x N /b;0) 
-p(^E|-Em| (53) 

Here b{t) is the scale factor satisfying the following equa- 
tion and initial conditions: 

b + uj 2 (t)b = ujl/b 3 , 6(0) = 1, 6(0) = 0; (54) 

Iq = \Jhj m.ujQ is the oscillator length. We assume that 
for t < to the frequency of the trap was fixed at u> = 
wo and Ej are the single particle energies corresponding 
to this frequency. The time parameter r(t) is defined 
according to r(i) = J Q dt' /b 2 (t'). 

Within this approach one can evaluate correlation 
functions as well as average energy for time-dependent 
trap. Here we consider the time-depending process of 
"switching off" the trap potential according to the (rela- 
tively general) law 

w 2 (t) =ojl(l-5t) r (55) 

where r is an arbitrary power < r < oo and 5 is the 
rate of the process. Note that r = 1 corresponds to 
the linear ramp considered in the main paper and r — > 
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oo with 8r kept constant corresponds to the exponential 
decrease of lo 2 with time. We will analyze the energy of 
the system at t = 1/8, i.e. when the trapping frequency 
vanishes: £(S,r). Note that in the limit 8 — > we must 
have £ (0, r) = because the gas can occupy the infinite 
volume. 

Using the scaling approach we described above one can 
relate the non-equilibrium energy density to the equilib- 
rium one via 



£(t,8 ) r) = ~£ CI 



8t) r b 2 (t)+b- 2 (t) + [b(t)} 2 (56) 



where £o is the initial energy of the Tonks gas in the trap 
and b(t) is a solution of the Eq. (|54p with Lo 2 (t) given in 
Eq. l[5"5"j). For simplicity we take w§ = 1. The equation 
for &(£) with (D 2 (i) = (1 — 8t) r can be solved analytically 
using the Ermakov approach (33|: the solution of the 
equation for b(t) is given by 



co (t)=exp(-<?Q 




0.000 0.002 0.004 0.006 0.008 0.010 

Ramp Rate 8 



b(t) 



(57) 



where xi^(t) are solutions of the linear Hill-type equa- 
tions 

x 2 +Cj 2 {t)x 2 {t) = 0, ar 2 (0)=0 ) i 2 (0) = 1 (59) 

These equations can be further solved in terms of the 
Bessel functions. However, their explicit form is rather 
cumbersome and we rather not show them here. We 
found that the the energy the rate 8 satisfies the following 
scaling 



FIG. 4: Energy of the Tonks gas after release from a trap as a 
function of the ramping parameter. The trapping frequency 
is changing according to Eq. ([55} with Wo = 1. Shown curves 
correspond to r = 1,2, oo, where the latter corresponds to 
the exponential decrease of lo 2 with time. All three curves are 
perfectly fitted by the power law dependence S oc S r ^ T+2 \ i.e. 
+ uj 2 (t)xi(t) = 0, £i(0) = 1, ii(0) = 0, (58) |5| 1/a , S 1/2 , and 5 for the three different curves respectively. 



£{8,r)^u; C{r) 



_8_ 



(60) 



where C(r) is a number of the order of unity. Note that 
at r — > oo the exponent in the power of \S\ saturates at 
unity, i.e. the system remains in the regime B). In Fig.fi] 
we illustrate the scaling (|60|) for several values of r. 

Thus, according to the classification scheme of the 
main text [l[ , the TG gas adiabatically released from the 
harmonic trap follows regime B). We note that in the case 
when the harmonic trap is not switched off completely or 
if one considers the process where the trapping frequency 
increases in time the energy dependence on cS becomes 
quadratic. So according to our classification the system 
is then in the analytic A) regime. In Fig. Owe show the 
corresponding dependence of A£(8) = £(8) — £(0) on 8 
at t — > oo for the trapping frequency increasing in time: 
u> 2 {t) = 1 +tanh(<5i). Note that not only the dependence 
is quadratic the heating is almost negligible compared to 
what one gets in the B) regime (see Fig. [4]). 

It turns out that the analysis above can be immedi- 
ately generalized to another very well known Calogero- 
Sutherland model describing one-dimensional Fermions 
interacting via 1/x 2 potential This model is a rear 
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FIG. 5: Excess energy of the Tonks gas after ramping on the 
trap frequency according to uj 2 (t) — l+tanh(<5£) as a function 
of S. The dependence perfectly agrees with the quadratic law 
AS oc S 2 , i.e. this process belongs to the analytic (A) regime 
according to our classification. 



example of solvable models which describe long-range in- 
teracting quantum systems. As an effective model it re- 
ceived a number of applications in various fields includ- 
ing quantum Hall effect, random matrix theory, etc. In 
a time-depending harmonic potential the corresponding 
time-dependent Schrodinger equation describing assumes 
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the form 
.5* 



1 dt ~ \2^dx 2 



N 

E 

j>i= 



A(A-l) 



W 2 (i) 



N 



E^ 2 r ( 61 ) 

3=1 



Using the scaling ansatz similar to the one employed for 
the Tonks gas, the particular solution of this equation, 
which corresponds to the ground state at t = 0, can be 
written as 1341 



*({*;},*) = 



1 



&iV/ 2 



JV 
J>i=l 



(62) 



where the function |6| satisfies Eq. ([54")) with the same ini- 
tial conditions. We again assumed w 2 (0) = 1 for simplic- 
ity. Note that there is close analogy between the scaling 
approach for the Tonks and Calogero-Sutherland models 
coming from the fact that the interaction energy in the 
latter scales in the same way with x as the kinetic energy. 
Similarly the " equipartition theorem" is satisfied for the 
model (|6T|) in equilibrium as well as for the TG gas: the 
half of the total energy comes from the harmonic well. 
This immediately implies that the scaling of the energy 
with 5 for the adiabatic turning off the potential is iden- 
tical for the two models. Thus we conclude that for the 
dependence J 1 = (1 — St) r the residual energy at t = 1/5 
again scales as |(5| r /( r + 2 ) and the model belongs to the 
class B) according to our classification. It is interest- 
ing that this conclusion (which is valid only for diagonal 
correlations, e.g. for the energy) does not involve the 
dependence on the parameter A. 



VIII. DICKE MODEL 

In this section we briefly consider a Dicke model, de- 
scribed by the following Hamiltonian: 



H 



Dicke 



= /x(t)ata + -JL [a f Sr + <*S+] , (63) 



JV 



i=l 



where a and a* are the bosonic fields and S + and S~ are 
the bosonic spin-1/2 razing an lowering operators. This 
Hamiltonian appears in many contexts of condensed mat- 
ter, atomic and optical physics. For example, this model 
represents interaction of N two level system with a pho- 
ton cavity mode. This model also describes Feschbach 
resonance of interacting fermions in the case of a broad 
resonance [351 ] . 

As a particular example we assume that fi(t) = —25t 
linearly changes in time and that initially (at t — > — oo) 



all spins are aligned along the z direction and the bosonic 
mode is empty. This setup is identical to that considered 
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FIG. 6: Bosonic occupation number for the Dicke model as 
a function of driving parameter 8 for different values of total 
number of particles N. The data is read off after sufficiently 
large evolution time. 



recently by A. Altland and V. Gurarie j35|. In the limit 
6^0 one expects that the system will follow the ground 
state and all spins will flip so that the population of the 
bosonic mode n = (a' a) at t — > oo is exactly N. In 
Ref. (35j it was indicated, however, that one approaches 
this limit in a nontrivial way. 



In Fig. ^ we plot the numerically found dependence 
1 — n/N as a function of the parameter 6 for various 
values of N. We were using a semiclassical approach 
similar to the one described in Appendix |Vj Let us first 
note that the dependence of n on 5 is linear in agreement 
with the regimes B) or C) in our classification scheme. 
Second we observe that there is no adiabatic limit for 
N — > 00 suggesting that in fact the dynamics belongs to 
the regime (C). We will present a more detailed analysis 
of the slow dynamics of the Dicke model in a separate 
publication [361 ]. 
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